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Using binding free energy to guide ligand design 
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The molecular distributions obtained from canonical Monte Carlo simulations can be used to find an approximate 
interaction energy. This serves as the basis of a method for estimating the binding free energy for a ligand to a 
protein which enables the free energy to be used to direct the design of Ugands which bind to a protein with high 
affinity. 
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I. INTRODUCTION 

The goal of structure-based drug design is to create a ligand 
which binds with high affinity to a protein target. An exciting 
prospect is the ability to carry out this design process com- 
putationally and thereby obtain a series of drug leads which 
are potent and which can be subsequently optimized for drug- 
like properties. In an earlier paper, we presented the worm- 
hole method 1 1] which, subject to some limitations, allows the 
binding affinity of a given drug ligand to a protein to be com- 
puted. For the purposes of structure-based drug design, we 
might imagine using the wormhole method to screen a large 
number of known molecules against the protein. This suffers 
from two serious drawbacks: only a tiny fraction of feasible 
drug-like chemicals can be assessed in this way; and it is not 
initially known where or how the compounds are likely to bind 
to the protein. 

One strategy for alleviating these problems is to use a 
fragment-based approach |2]. Here fragments are small (usu- 
ally rigid) organic molecules. By a judicious choice of frag- 
ments, a large and diverse set of drug-like molecules may be 
built in silico by forming bonds between them. Because the 
number of fragments is small, O(IOO), and because they are 
relatively simple, we can compute maps of where the frag- 
ments bind to the protein. These data then serve as the build- 
ing blocks to create larger drug-like molecules. The process of 
constructing these molecules, of necessity, provides the bind- 
ing mode to the protein. The key is, of course, to build the 
large molecules in such a way as to optimize the binding affin- 
ity. 

We face two challenges here. Given a partially built candi- 
date molecule, can we quickly assess how a particular frag- 
ment can be added? Having grown the molecule with the 
addition of a fragment, how can we rapidly compute the re- 
sulting binding affinity to evaluate whether the new molecule 
is acceptable? 
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In this paper we describe possible solutions to both these 
problems. Let us start by giving an overall description of the 
method. Our computational model consists of a protein which 
is either kept rigid or is allowed to have a small number of 
degrees of freedom interacting with a ligand. We assume that 
the ligand is "made" by connecting several simpler organic 
fragments together and, for simplicity, we take the fragments 
to be rigid and assume that they are joined by rotatable bonds; 
however, the method is easily generalized to remove these re- 
strictions. This system is described by a conventional force 
field such as Amber |3] with the effects of the solvent cap- 
tured by an implicit solvent model such as GB/SA |4; 5]. In 
this model, we can limit the number of degrees of freedom of 
the system to a manageable number, 0(10). For the purposes 
of this discussion, we assume that we have identified a binding 
site on the protein. Our goal is to design a set of ligands (cre- 
ated from the fragments) which bind to the protein with high 
affinity. It would be possible to build additional criteria into 
the design process, e.g., synthesizability, solubility, etc.; how- 
ever, these considerations are beyond the scope of this paper 
Our standard for success is that the approximate techniques 
we develop lead to ligands with high affinity as predicted by 
the full force field outlined above. For this purpose it is conve- 
nient to regard the full computational model as "exact". The 
degree of agreement with experimental data, while crucial, en- 
tails validation of the force field which, again, is beyond the 
scope of this work. 

We begin by performing wormhole Monte Carlo simula- 
tions (3| of each of the various fragments binding to the pro- 
tein. These calculations give the binding affinity of each 
fragment to the protein and equilibrium distributions of the 
fragment-protein system. We next fit an analytic function, a 
Gaussian mixture, to these molecular distributions. The fits 
for two fragments are then used for two purposes: to find fea- 
sible ways to form a bond between the fragments, creating a 
larger ligand and to give an approximate interaction energy 
for the newly created ligand with the protein which, in turn, 
allows for rapidly computing its binding affinity via the worm- 
hole method. 

The first part of this paper describes techniques for fitting 
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a Gaussian mixture to a distribution of molecular configura- 
tions. We adopt the well-known EM method [6J for this pur- 
pose; however, we have to adapt the method to deal with two 
peculiarities of molecular distributions: firstly, in the presence 
of constraints, the distributions lie on a sub-manifold of Carte- 
sian space; secondly, for symmetric molecules, we can make 
a better fit by respecting the symmetry. 

Two applications of Gaussian mixtures for molecular distri- 
butions are described next. They may be used to define suit- 
able portals for the wormhole method; this provides a more 
robust method than the use of ellipsoidal portals given in |l|]. 
They may also be used to provide an approximation for the 
energy of a molecular system. 

Finally, we describe how these tools may be combined to 
compute an approximate binding affinity which allows the 
binding free energy to be used to direct the design of ligands. 

II. GAUSSIAN MIXTURES 

There is often an interest in fitting some observed data with 
a "model", an analytic function which approximates the data. 
One important category of data is the set of configurations 
of a molecular system given, for example, by the results of 
a Monte Carlo simulation. An analytic fit then provides an 
approximate but compact representation of the observed data. 
Because the samples from a canonical ensemble Monte Carlo 
simulation are drawn from a distribution which is proportional 
to cxp(— /3i?(x)), where -E(x) is the energy of the system in 
configuration x, /3 = l/{kT), k is the Boltzmann constant, 
and T is the temperature, the analytic fit can also be used to 
give an approximate expression for the energy of a molecular 
configuration. 

An important class of models is the mixture of Gaussians 
and the EM (expectation-maximization) algorithm |6] is fre- 
quently used to optimize this model based on the maximum 
likelihood. We begin by reviewing an iteration of the standard 
EM algorithm including the straightforward extension of al- 
lowing the samples to have a statistical weight. Assume that 
our data is 

[xi,X2,X3, . . . ,X„], 

where Xi is a point in R'' and that associated with each of the 
samples x^ is a scalar weight Wi. This weight might arise from 
coalescing consecutive identical samples from a Monte Carlo 
simulation (because of a run of rejected moves) or because 
the Monte Carlo sampling is carried out with a non-physical 
energy E* in which case we have Wi — exp[—P{E{xi) — 

Let the current fit be 
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ajG(x;yj,Cj), 



where Y^^=o^ '^j — 1 ^^'^ ^(x;yi,Cj) is a d-dimensional 
Gaussian with unit volume and mean, yj, and covariance, Cj. 
The goal is to find the set {aj,yj, Cj} which maximizes the 



log-likelihood 

L = (ln/(xj); Wi)i, 
where (•; •)i denotes the weighted arithmetic mean. 
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An EM iteration proceeds as follows: 
ajG(xj; yj,Cj) 
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With this procedure, the log-likelihood is guaranteed to con- 
verge to a local maximum ja]. The quantity gij gives the re- 
sponsibility of the Gaussian component G(-; yj, Cj) for x^. 



A. Non-Cartesian space 

When simulating complex molecules it is important to re- 
duce the dimensionality of configuration space by imposing, 
for example, constraints on bond lengths and bond angles. For 
example, when simulating biphenyl (two benzene rings linked 
by a single bond), the energetics of the molecule allows us to 
treat it as two rigid rings connected by a bond that permits 
only torsional movement. The complete configuration of the 
molecule is then given by the position and orientation of one 
of the rings together with the torsion angle of the connecting 
bond. 

We will represent torsion angles as a point on the circle 
Orientations are conveniently represented as unit quaternions 
0]; however, because q and —q represent the same orienta- 
tion, orientations are defined as a pair of opposite points on 
S'^. In Mardia and Jupp |8], a distinction is made between a 
directed line through the origin, a direction (which can rep- 
resent the torsion angle), and an undirected line through the 
origin, an axis (which can represent orientations of general 
molecules). The orientation of a diatomic molecule, for ex- 
ample CO, would be given by a unit vector, i.e., a direction 
on The full configuration x is then a mixture of Cartesian 
coordinates and "angle-like" coordinates in Our strategy 
for applying Gaussian fits to points in this mixed topology is 
to replace eq. ( l3ct by 

yf- = ((x,;g,,i«,))„ 

where {{■;■)){ is the appropriate weighted "physical" mean of 
Xi. For the angle-like coordinates, we find the mean by em- 
bedding §' in M'+i. The mean direction is given by the di- 
rection of the weighted sum of the unit vectors while the 
mean axis is given by the axis about which the moment of 
inertia of the weighted axes is minimum f7;'^. 
Similarly, we replace eq. i3d\ by 
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where d(x, y) is a displacement in from y to x. The op- 
eration of d is to map configurations into a local Cartesian 
space centered at y. In order to make such a mapping for 
the angle-like coordinates, we project the sphere (in the case 
of directions) or hemisphere (in the case of axes) onto a ball 
in M' using a generalization of the Lambert azimuthal equal- 
area projection |1] with the pole of the projection given by 
the mean. It is important that the projection preserve area so 
that its Jacobian is constant; in this way, integrals in the pro- 
jected space are the same (up to a multiplicative constant) as 
integrals in the original space. 

In carrying out this extension of Gaussians to non-Cartesian 
geometries, we have lost an important property of the Gaus- 
sian. In R*^, if we fit a single Gaussian to arbitrary data, then 
the log-likelihood is maximized by choosing the mean and co- 
variance of the Gaussian equal to those of the data. We are not 
aware of a generalization of the Gaussian which preserves this 
property for our more complex geometries. However, the pre- 
scription given above presumably nearly preserves this prop- 
erty provided that the covariances of the individual compo- 
nents in the mixture are sufficiently small that the Gaussians 
do not "wrap around" §' to any great degree. We will address 
this issue later. 

In the discussion above, we have implicitly assumed that 
a uniform measure on is the natural metric for angle-like 
coordinates. This is the case for the orientation of a mole- 
cule and for torsion angles. However, the situation is more 
complex for molecules whose bond angles can vary (for ex- 
ample to treat the common conformations of cyclohexane). A 
full treatment of such cases is beyond the scope of this paper. 
However, the strategy would be the same as given here: de- 
termine a suitable mean and then map the samples to a locally 
Cartesian space centered at the mean in such a way that con- 
figuration space integrals can be expressed in the transformed 
space with a constant Jacobian. 

B. Incorporation of symmetries 

The use of symmetries allows the simplification of many 
problems. In describing molecular configurations, we en- 
counter both discrete and continuous symmetries. Examples 
of the latter are translational and orientational invariance when 
simulating a solute molecule in a large volume of solvent, ro- 
tational invariance about the axis in a diatomic molecule, etc. 
Such symmetries are best treated by expressing the molecu- 
lar configurations in a lower dimensional space thereby ig- 
noring the symmetry coordinates. Thus the "orientation" of a 
diatomic molecule can be expressed as a direction on rather 
than as an axis on S"^. 

Let us describe some typical discrete symmetries that arise 
in molecular systems. A molecule of methane, CH4, may 
be oriented in 12 different ways (the order of the tetrahedral 
group T) that leave like atoms in the same positions. We do 
not treat the reflection symmetry of methane as an additional 
symmetry because such inversions do not occur under normal 
conditions. 

A more complex example is biphenyl. When bound to a 



protein, this has 8 symmetries made up of combinations of 
180° rotations of the benzene rings about the connecting bond 
and an interchange of the two rings. When biphenyl is placed 
in any of its 8 symmetric positions the resulting system en- 
ergy and hence the equilibrium distribution is the same. If, 
on the other hand, the biphenyl is free in solution, then we 
remove the continuous symmetries by fixing the position and 
orientation of one of the rings. There are then 4 symmetries 
given by ijj ±ip and ip ^ ±'ip + tt where ip is the torsion 
angle. These correspond to rotating the free ring by 180° and 
changing the sign of the torsion angle. The latter operation 
places the biphenyl into its mirror symmetric conformation 
(but not by inverting the molecule). This symmetry is nor- 
mally excluded when biphenyl is bound to a protein, because 
the protein binding pocket will not exhibit the same symmetry 
(because proteins are chiral). 

We shall suppose that the system has a fc-order symmetry, 
which can be described by a symmetry operator S{-, I) where 
<l < k and 

X 1-^ S'(x, I) 

maps the configuration into one of the k symmetric con- 
figurations. We can compose symmetry operations with 
S{S{xJ')J) = S{x,l © /'). Clearly © defines a group of 
order k. We will take identity element to be and define the 
inverse of Z to be ^ (thus, Z © Z = 0). 

In fitting a Gaussian mixture to data, we can use S both to 
symmetrize the samples and to symmetrize the fit. However, 
by using the properties of S the computational complexity in- 
creases only by k (instead of k^). We begin by symmetrizing 
the fit, 

m— 1 k — 1 

/(x)= ^a,-^G(5(x,0;y„C,). 

j=o 1=0 

For simplicity, we apply the symmetry operation through the 
configuration argument of G rather than via yj or Cj . From 
this definition, it is easy to show that 

/(5(x,0) = /(x). 

(This follows from the group properties of ©.) In forming the 
responsibility matrix, we start by computing the responsibility 
of the component G{S{-, l');yj, Cj) for the symmetrized data 
point S{xi,l), 

k /(5(x„0) 

1 ajG'(5'(x,,Z^©/);y^-,Cj-) 

k /(xO 

= 9ij(l'@l) 

where 

_ 1 ajG'(5(x,,0;yj,Cj) 
~ k /(xO 
We can now update the components using 

new / \ 

Qf-^ = (d(5(x„/),yf-)d(5(x„0,yr")^;g.,i«;»)»,i, 
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Here in forming 

I E [0,k). 



and 



), we sum over i G (0, n] and 



The partial EM update then consists of updating the respon- 
sibiHties for the new component. 



C. Extension of the greedy algorithm 

In the foregoing, we have supposed that the number of com- 
ponents in the fit is known. In general, this is not the case and 
various algorithms have been proposed to grow the number of 
components in such a way that a fit close to the global max- 
imum for the log-likelihood is tracked. Here we adapt the 
greedy EM algorithm |9] for adding components so that sym- 
metries can be included. We determine the optimal number 
of components by minimizing a cost function involving the 
minimal description length llOl §7.4.2], 



(4) 



The term in brackets gives the number of free parameters in 
an m-component fit and p, which is normally unity, is a pa- 
rameter that can be adjusted to penalize the addition of more 
components. 

Let us review the greedy algorithm f^. After the EM al- 
gorithm has converged for an m-component fit, we attempt 
to add a new component (with index m) as follows. Initially, 
each data point is assigned to the component j for which 
gij is maximum. In this way the data is partitioned into ni 
sets Aj . We make several splits of each component j by se- 
lecting two random samples from Aj and partitioning Aj into 
two subsets based on closeness to the two random samples. 
A tentative new component is added with a,„ = aj/2 and 
mean and covariance given by one of the two subsets. The 
resulting tentative fit undergoes partial EM iterations where 
ctm, ym, and Cm are adjusted and the aj for < j < m 
are merely scaled by 1 — a„i (with the corresponding y-j and 
Cj held fixed). This procedure is repeated several times for 
each of the m components and the fit with the maximum log- 
likelihood (following the partial EM updates) is selected as 
the (to + l)-component fit which is then subjected to full EM 
updates. 

When treating weighted samples, we modify the procedure 
above by selecting the two components from Aj with proba- 
bilities proportional to their weights. Because we do this sev- 
eral times, we use the Walker algorithm 1 1 1 1 to make these 
selections. 

In order to include symmetries, we generalize Aj above to 
Aji' which contains those 5(xi, I) for which gujii — gij(VQi) 
is maximum. It is only necessary to consider splitting the 
TO unsymmetrized components of the existing fit; thus we 
only need to determine AjQ. We can do this by assigning 
an unsymmetrized sample to a symmetrized component by 
finding the j and V which maximizes gijii, and then adding 
S{-K,,l = l')toAjo. 

As before, we partition each AjQ into two subsets by pick- 
ing two random samples from (according to the sample 
weights) and using the distance as defined by d as the close- 
ness metric. For each subset we use an initial = 1"^ and 
(ym, Cm) computed from the data in the subset. 



amG{S{:x.i, l);ym, Cm) 



(1 — am)kf{-x.i) + arn G(S'(xi, Z'); 

ym ; Cn 



where we need to evaluate gimi for I G [0, k) and for i e AjQ, 
i.e., for all i, for which S{xi, I) S Ajo for some /. The update 
of the TOth component is then 



now / \ 

"m = {9imi;Wi}i',l, 



- (d(5(x„0,y:^r)d(5(x„/),y,"r)^;ff.m;w'.).-,;, 

where the subscript i* indicates that the sums over i should 
include only i £ Ajo. 

D. Loose ends 

Finding the 1 -component fit with non-symmetric data is a 
simple matter of computing the mean and covariance of the 
data. However, if we are performing a symmetric fit, we need 
to apply EM iterations to obtain a converged one-component 
fit. To determine a starting point for these EM iterations we 
pick a random "central" sample, and transform the other sam- 
ples using the symmetry operator so that they are as close as 
possible to the selected sample. The resulting n symmetry- 
transformed samples are used to define a tentative (yo,Co). 
This procedure is repeated several times with different cen- 
tral samples and the (yg, Cq) which yields the maximum log- 
likelihood is used as the initial guess for the first component. 

The EM algorithm can fail with poorly conditioned sam- 
ples. For example, one Gaussian component might converge 
to a group of samples which are in a lower dimensional space. 
We avoid this problem by placing a lower limit on the max- 
imum eigenvalue of the covariance matrix and by placing a 
lower limit on the ratio of the minimum to maximum eigen- 
values. This makes the EM algorithm more robust possibly at 
the cost of requiring more components to maximize the log- 
likelihood 

In Monte Carlo applications, we may wish to avoid Gaus- 
sian components where any of the angle-like coordinates 
"wrap around". The presence of such wrapping may destroy 
detailed balance because a transition to a wrapped sample 
drawn from such a Gaussian is not balanced by a reverse pro- 
cess. We can limit the effect of the wrapping by checking 
those diagonal elements of the covariance matrix correspond- 
ing to the angle-like coordinates. If these are so large that 
wrapping occurs within 3 standard deviations of the mean, 
for example, then we can scale the corresponding rows and 
columns of the covariance matrix appropriately so that wrap- 
ping is limited to the small fraction of samples beyond 3 stan- 
dard deviations. Here again, the algorithm can adjust to this 
constraint with additional Gaussian components. In the next 
section, we will show how detailed balance can be maintained 
exactly for wormhole moves even in the face of wrapped 
angle-like coordinates. 
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The result of a canonical Monte Carlo simulation is a set 
of configurations Xj drawn from the underlying Boltzmann 
distribution proportional to exp(— /3i?(x)), together with the 
corresponding energies E{xi). In the foregoing discussion, 
we make the fit to the configurations, essentially ignoring the 
energies. This is an appropriate use of the data from a Monte 
Carlo simulation where the sample configurations constitute 
the "primary" data. One application where we could make 
use of the energies is when making fits to several indepen- 
dent Monte Carlo runs of the same system. In this case, we 
can adjust the overall weight of each independent run so that 
the difference of f3{E{xi)) and (ln/(x,;)} is approximately 
the same across the runs. (Here, (•) denotes a average over a 
single run.) This adjustment is important when the individual 
runs are not sufficiently long to sample configuration space 
fully. 



III. APPLICATIONS OF GAUSSIAN MIXTURES 

We use the procedure for fitting molecular distributions 
with a Gaussian mixture in two ways. The first is asa method 
of defining the portals for wormhole Monte Carlo 1 1 ] . In this 
case we are fitting the data from several independent Monte 
Carlo runs and Gaussian mixtures then offer a robust way of 
"clumping" the data with each component of the mixture then 
providing a portal for the wormhole method. The other appli- 
cation provides an approximation to the energy of the system. 
Here we are more concerned about the accuracy of the fit, and 
we also need to establish that the arbitrary constant that con- 
nects the energy to the logarithm of the fit drops out when 
forming physically relevant quantities. 



A. Portals for wormhole Monte Carlo 

The original description of wormhole Monte Carlo 1 1 ] was 
specialized to the treatment of molecular dissociation, 

A + B ^ AB, 

where we sought the equilibration between the bound and un- 
bound states of molecules A and B. This procedure can be 
generalized to deal with other types of interaction, e.g., molec- 
ular exchange. 



AB + C^ A + BC, 



protonation. 



A + H^ 



AH+, 



or tautomerization. 



ABE ^ HAB. 

(In practice, the free proton in the second case would be han- 
dled by an implicit solvent held at constant pH.) We therefore 
consider the equilibrium of A "systems" indexed by A. Each 
of the systems is made of ^\ independent molecular "com- 
plexes" indexed by (p and each complex is made up of 1 or 



more interacting molecules. (Thus with molecular dissocia- 
tion we have A = 2. The unbound system A = consists of 
$0 = 2 independent complexes, each consisting of a single 
molecule, A or B, while the bound system A = 1 consists of 
$1 = 1 complex, AB.) 

If the configuration of complex (j) in system A is k\^, then 
the full phase space is given by T = {A; x^o, x^i, . . . , }. 
Here, we have added a set of "ignorable" coordinates, x^o; the 
energy of the system, and hence the equilibrium distribution 
function, is strictly independent of these coordinates. For ex- 
ample, in simulating a molecule in a solvent bath, x^o would 
include the position and orientation of the molecule; or, when 
a molecule is deprotonated, it would include the coordinates 
of the "missing" proton. Inclusion of these ignorable coor- 
dinates is dictated by the requirement that T span the same 
phase space volume for each A. In practice, we do not keep 
track of x^o, because the integral over this coordinate is trivial 
(the integrand is constant!) and we write 

dxAo = vxQ. 

Wormhole Monte Carlo moves allow the state to switch be- 
tween different systems preserving detailed balance. This al- 
lows the determination of the ratios, 

Wo:Wi:W2:..., 

where Wfj, is the statistical weight of system fi. 



exp(-/3F^)= J Sx^eM-PVdT, 



and Ffj_ is its free energy. In particular, in the case of protein- 
ligand binding, the dissociation constant is given by 

where Vq is the system volume. 

In this more general framework, the wormhole move [IJ is 
defined as follows. We define a set of "portal functions," w, 
w', w", . . . , on T, with properties 

< w(T) < 1/v < oo, 
/dTw(T) = 1, 

where u is a representative phase-space volume of the portal 
function. A wormhole move consists of the following steps: 
select a pair of portals {w, w') with probability Pww'', reject 
the move with probability 1 — vw{T), where T is the current 
state; otherwise, with probability ^^(T), pick a configuration 
T' with probability w'(T'); and accept the move to T' with 
probability 



Pw(T,T')=min 1 



p^,^ exp{-PE*{V))v' 
' Pww' exp(-/3i?*(T)) V 



(5) 



where £'*(T) « E{T) is the "sampling" energy of configu- 
ration T, which is used also for the conventional Monte Carlo 
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moves (within a system). We term w and w' the source and 
destination portals, respectively. The test involving w{T) de- 
termines whether the current configuration is "in" the source 
portal — note, however, that this test is "fuzzy". If the test suc- 
ceeds, a move is attempted to the destination portal, and the 
move is accepted according to a standard Boltzmann factor 
modified by the ratio of the portal volumes. In the limit of a 
long Markov chain, we then have 

-> C{Sx, exp{-P[E{T) - E*{T)])), 

where C is independent of ^ and (•) is the average over the 
Markov chain. 

Although the choice of portal functions is arbitrary, the 
method is only effective if vw{T) is sufficiently large to al- 
low wormhole moves. For simplicity, we restrict each portal 
function to a particular system A = /i. Because the complexes 
making up a system are independent it is natural to consider 
^(T) as product of density functions for each complex. Thus 
the typical portal function is 

1 

where 

J Wfj,^{Xf,^) dx^,0 = 1, 

and the factor l/Vf^o arises from an implicit constant den- 
sity in Xpo- For a particular complex (f) in system A we need 
to determine a set of portal functions ^^^^(xp^), u>^^(xp0), 
w^^(x^0), . . ., which reflect the probable configurations for 
this complex. We obtain these portal functions using the re- 
sults of several conventional canonical Monte Carlo runs on 
the complex. We make a Gaussian fit to the resulting sets of 
configurations. If the fit contains m components, then we ob- 
tain m portal functions, indexed by j, for this complex each 
of which is a symmetrized Gaussian of the form 

where kx^ is the symmetry order for the complex, 5*^0 is the 
corresponding symmetry operator, etc. We take the "volume" 
of this portal function to be 

vxc/ij — kx4>IG{yx4>j\y\ci>jT^X4>'j)- 

We assume that the choice of source and destination portals 
is independent so that the portal probability Pww' can be fac- 
tored into probabilities for w and w'; furthermore, we assume 
that these probabilities may in turn be factored into choices 
for the source and destination systems and for the portals for 
the respective complexes for each system. In this case, the 
wormhole moves can implemented as follows. Pick a source 
portal system /i; if A ^ /i, the move fails; otherwise consider 
each complex in the system /i in turn; for complex 0, pick a 



random portal function j and pick a random symmetry index 
I; with probability 

1 ^ 'yp0jG(S'p0(x^0, Z); yfj,ipj, C^^j)/^^^, 

reject the move; if none of these tests cause the move to 
be rejected, the "in" test succeeds and we proceed with 
choosing the destination portal by picking the destination 
system /i', picking a portal function j' and a symmetry I' 
for each complex 0', and setting the configuration for the 
complex to Sfj_'^'{x^'^',l') where x^/^/ is selected from 
G(x^/0/ ; Yu^'f , )■ In evaluating the acceptance prob- 

ability, eq. (|5), we express Pww' as the product of the individ- 
ual probabilities (of selecting source and destination systems 
and of selecting particular portals for the source and destina- 
tion complexes). Similarly the volume of the portal is given 
by the product of v^o and the volumes of the portal functions 
for the separate complexes. 

In this formulation, the symmetry of a complex is incorpo- 
rated into the portal function, eq. (|6j. However the test for 
being in the portal and the operation of selecting a configura- 
tion from it are decomposed into picking a random symmetry 
(with equal probabilities) followed by a test or selection on a 
unsymmetrized Gaussian. 

If, when sampling from the destination portal, any of the 
angle-like coordinates are wrapped around, then we immedi- 
ately reject the whole move. This is necessary in order to 
maintain detailed balance, because the test on the source por- 
tal never involves wrapped coordinates. This effectively re- 
places the Gaussians in the definition of the portal functions 
by clipped versions which evaluate to zero for wrapped coor- 
dinates. 

There is a great deal of flexibility in the choice of portal 
probabilities offered by the scheme outlined here. Because 
the test of being in the source portal is typically very inex- 
pensive, it is desirable to arrange that the source portal prob- 
abilities are roughly equal. In addition, we usually adjust the 
ratio of conventional to wormhole moves so that, on average, 
each configuration is tested against all the portals for every 
attempted conventional move. On the other hand, the proba- 
bilities for the destination portals would usually be adjusted to 
reflect the statistical weight of the portal. 

Let us turn to the details of making the Gaussian fits to 
define the wormhole portals. Because the individual Monte 
Carlo runs performed for each complex are independent, it 
is natural to consider scaling the overall weight to the results 
from each run in order to match the energy samples. In prac- 
tice, this procedure results in rather poor fits with too many 
components. In this application, Gaussian fitting may viewed 
merely as a robust clumping technique and for this purpose if 
suffices to attach the same weight to all the samples. For the 
same reason, we increase p to 5 in eq. (|4} so that a smaller 
number of components is used to make the fit. 

Gaussian portals offer advantages over the use of ellipsoidal 
portals proposed in llj. With a given number of components, 
the EM method does a "global" optimization and is thus likely 
to obtain a better fit and than the somewhat ad hoc scheme for 
choosing ellipsoids. Also the Gaussian fit to a configuration 
of independent complexes naturally factors into a product of 
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Gaussians for each complex. Thus Gaussian portals reflect 
the independence of complexes properly. Gaussian portals, 
combined with the mean energy for a portal (which can be 
readily estimated from the energies of the samples) also offer 
a rough a priori estimate of W\. This, in turn, allows us to 
adjust vxo to maximize the probability of transitions between 
systems and hence to reduce the error in the eventual estimate 
of W\. In the case of Ugand-protein binding, where we simu- 
late single ligand and protein molecules in a system of physi- 
cal volume Vq, this procedure entails adjusting Vq so that the 
fraction of time the molecules are associated is roughly |. 

Finally, we remark that when performing a conventional 
Monte Carlo move for a particular system, it is preferable to 
select randomly a single complex to move. This will result in 
a higher acceptance rate compared to attempting to move all 
the complexes simultaneously. 

B. Obtaining the energy from the fit 

The result of a wormhole Monte Carlo simulation is a set of 
configurations sampled from exp(— /3i?* (T)). If we fit an an- 
alytic function to these samples weighted by exp{—(3[E{T) — 
E*{T)]), then the fit can serve as a basis for approximating 
E{T). Here we detail how we can use Gaussian mixtures to 
carry out this fit and we show how to obtain approximations 
for the energies of the individual complexes and how the arbi- 
trary offsets for energies cancel whenever energy differences 
are computed using the approximate energies. 

We begin by making normalized fits, f\^, to all the com- 
plexes in all the systems. If the same complex appears in mul- 
tiple systems, the samples may be aggregated in order to per- 
mit a fit using more data. The energy of the system is taken to 
be the sums of the energies of the complexes, i.e., 

4,=1 

and similarly for E*{T). The sampled configurations for 
each complex are assigned weights of exp{—P[E{x\^) — 

The normalized fit for a particular system is then given as 
the product of the fits for the contributing complexes, multi- 
phed by 1 / vxo and these can be combined weighted by Wx to 
provide a fit in T space as 

Wx 1 

^" 0=1 

where W = J2x and we have 

-pE{r)^D + lnf{r), 

where D is an arbitrary adjustable constant. This provides an 
approximation of the energy of a system. In addition, we can 
approximate the energies of the individual complexes by 

-/3-E(xa0) ~ Dxcj, + ln/A0(xA0), 



where Dx^ are adjustable constants which satisfy 

^Dx^ = D + IniWx/W) - Invxo, 

</.=! 

for all A. 

Let us apply this to the case of molecular dissociation. The 
A = (resp. A = 1) system contains two complexes (resp. one 
complex) each of which is free to move within a system 
of 3-dimensional volume Vq; thus, we have vqo = (fVb)^ 
(resp. v\Q = (jVo), where a is the volume of orientation space. 
In this case, there is one constraint on the choice of energy off- 
sets for the fit energies, namely 

Dw = DQO + DQi-\Yi{Kd/a). 

As expected, this constraint does not involve Vq. It is also 
apparent from the form of this constraint, that differences in 
energy between the unbound and bound systems will be inde- 
pendent of the choice of offsets. 

iV. iUIOLECULAR DESiGN 

We now have the tools to tackle ligand design. We start the 
process by computing the binding affinity of the fragments. 
Fragment-based design works on the principle of building a 
complex molecule from simpler sub-components. We extend 
this idea by also computing the binding affinity of the larger 
molecule using data from the calculation of the binding affin- 
ity of the simpler molecules. Finally, we describe the process 
by which simple molecules can be combined. 

A. The single fragment binding affinity 

Our starting point is a protein target for which we know the 
structure and a set of simple organic fragments. The symme- 
tries of the fragment are determined. In the case of a rigid 
fragment, this consists of the set of 3-dimensional rotations 
which leave the molecule invariant. The energy of the system 
is computed using a conventional force field with an implicit 
solvent model as described in the introduction. The simula- 
tion is focused on a certain portion of the protein by adding a 
restraint energy which is zero if the fragment is within a re- 
gion of interest on the protein (e.g., within a binding site) and 
increases parabolically outside this region. It is possible to de- 
fine the restraint region to include a few solvent layers about 
the entire surface of the protein — ^but this obviously results 
in a longer simulation. Including a parabolic portion to the 
restraint potential allows the fragment distribution to fall off 
gradually and this allows the distribution to be fit with fewer 
Gaussian components than with a hard restraint. The binding 
affinity of a molecule will be only weakly dependent on the 
precise extent of the restraint region providing that it encom- 
passes the true binding site of the protein. 

We perform a wormhole calculation to find the binding 
affinity and to provide the distribution of fragments. In or- 
der to determine initial portals for this calculation we system- 
atically search for plausible binding modes by inserting the 
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fragment randomly into the restraint region with a random ori- 
entation and random conformations for the protein and frag- 
ment (if these molecules are flexible). This process is most 
efficiently carried out with a tailored restraint region (which 
prevents attempts to insert the ligand within the protein) fol- 
lowed by a quick steric check (where configurations are re- 
jected as soon as two clashing atoms are found). This can be 
followed by a crude energy minimization using the vacuum 
energy model. We can make an estimate of how many probes 
need to be made in order to explore the surface of the protein 
thoroughly and so to find all possible binding pockets. This 
estimate is based on the volume of the restraint region and 
the typical length and orientation scales for energy variation. 
A similar exercise is carried out for the unbound system — this 
merely consists in finding allowable conformations of the pro- 
tein and ligand. We then drop any of the bound configurations 
whose energy exceeds the minimum unbound energy. The re- 
sulting bound and unbound configurations are used as starting 
points for a set of conventional Monte Carlo runs with the 
full sampling energy. The initial portion of each run should 
be discarded and any bound run whose energy is stuck close 
to (or above) the unbound energy should be eliminated. The 
resulting data from these Monte Carlo runs is then used to 
determine wormhole portals using a Gaussian mixture and to 
estimate a starting value of the system volume Vq . 

In addition, we can add "catch-all" portals for the unbound 
and bound systems. For the unbound system, this will allow 
the molecules to assume arbitrary conformations (subject to 
whatever constraints are imposed by the molecular model). 
For the bound system, the molecules would be allowed to 
assume arbitrary conformations and in addition the ligands 
would be selected from the restraint region with an arbitrary 
orientation. These catch-all portals allow new binding modes 
to be discovered. 

The binding affinity is then calculated using the wormhole 
Monte Carlo method. During the course of this simulation, 
Vq is adjusted to maintain Wq ~ Wi and if Vq is increased 
(resp. decreased) we reduce Wi (resp. Wo) by the same factor 
We may also find that the ligand becomes trapped in a local 
energy well. Whenever this is detected (by the absence of 
successful wormhole moves), new Gaussian portals are found 
by rerunning the Gaussian fit adding recent configurations and 
restarting the binding affinity calculation. 

The process converges with a sufficiently long run with- 
out the need to add new portals and with sufficiently frequent 
wormhole moves between the bound and unbound systems. 
This process provides an estimate of the dissociation constant 
Kd for this fragment-protein interaction and a set of samples 
for the bound and unbound configurations. From this con- 
figurational data (weighted to reflect the difference between 
the full and sampling energies), we fit a Gaussian mixture to 
the bound and unbound distributions. This provides an ap- 
proximation to the energy of the protein and the ligand either 
unbound or as a bound complex. 



B. Approximate energy of combined molecules 

Let us consider the case where we have identified two pos- 
sible ligands AB and BC and we wish to combine these via 
the "overlap" portion B to form a ligand ABC. We might 
form a A^-fragment ligand with A and C being fragments and 
B being an {N — 2)-fragment overlap ligand. Alternatively, 
we might take B to be null and merely add a fragment C to 
a {N — 1) -fragment ligand A. (In either case, we form a 2- 
fragment ligand by taking A and C to be fragments and B to 
be null.) 

Because we are concerned here with the energies of differ- 
ent combinations of fragments, we adopt a notation for the 
energy where we do not explicitly specify the molecular con- 
figuration and where E{Z) is the energy for a single molecule 
Z and E{X, Y) is the energy for the two interacting molecules 
X and Y. We write 

£;(ABC) = £;(AB) + £:(BC) - £;(B) 

+ 6E{ABC), (7a) 
£;(ABC,P) = £;(AB,P) + £;(BC,P) - £;(B,P) 

+ 6E{ABC,P), (7b) 

where we expect 5E to be small if the energies are approx- 
imately additive. In this and subsequent equations, we un- 
derstand the configuration of the molecules to be consistent 
throughout the equation, e.g., A is in the same configuration in 
all the terms. If B is nufl, then we can write E{B,P) — E{P). 

We now make two approximations. We assume that the en- 
ergies involving the simpler molecules AB, BC, and B, are 
given by fits to configurations from prior binding affinity cal- 
culations and we assume that SE{ABC,P) and 5E{ABC) 
are small. If these assumptions hold, then eq. Q provides 
a method of computing the energy of the more complex mole- 
cule ABC very rapidly which, in turn, allows its binding affin- 
ity to the protein to be determined quickly. We expect the 
neglect of SE{ABC, P) and 6E{ABC) to be most easily jus- 
tified when the overlap portion B is as large as possible; i.e., 
when two {N — 1) -fragment ligands are combined to form 
an iV-fragment ligand. In carrying out this calculation, the 
arbitrary constants that enter when converting the fits to ener- 
gies cancel when considering energy differences between the 
bound and unbound systems. 

In our initial implementation, we build ligands by adding 
fragments one at a time with no overlap portion (i.e., B is 
null). Furthermore the approximate expressions for the energy 
are applied recursively so that the energy of a iV-fragment 
ligand interacting with the protein is found by summing each 
of the individual fragment-protein energies. 

The accuracy can be improved as follows: Assume that we 
have computed the binding affinity of the best ligands with up 
to — 1 fragments using the full energy. Fits to the distribu- 
tions of these molecules provide the corresponding approxi- 
mate energies which can be used to compute the approximate 
energy, using eq. Q, for A^-fragment ligands either by adding 
a single fragment or, preferably, by using an ( A^ — 2) -fragment 
overlap. This allows us to compute approximate values for the 
binding affinity of the A^-fragment ligands. The binding affin- 
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ity of those ligands with the best approximate binding affini- 
ties can then be recomputed using the full energy. Because 
this latter binding affinity calculation is carried out following 
the similar calculation with the approximate energy, we can 
use Gaussian fits to the samples from the approximate calcu- 
lation to define the initial portals for the wormhole method 
with the full energy. This procedure can be repeated to create 
ligands with an arbitrary number of fragments. 

In order to compute the full energy of the enlarged mole- 
cule, we may have to determine the Amber atom types afresh, 
for example, using the GAFF rules |12; 13]. In addition, 
we need to determine the partial charges for the new ligand. 
One simple prescription is as follows: when two fragments, 
A and C are combined to form a 2-fragment molecule AC, 
the charge on the hydrogen removed from A is donated to the 
derivatized atom on C and vice versa. This rule, which main- 
tains charge neutrality, can be readily generalized for larger 
ligands. More realistic charge models could be employed, if 
necessary. For example, VC/2003 1 14] allows partial charges 
to be computed without the need for a quantum calculation, 
while AMl-BCC \ l3,^ gives the charge on the basis of an 
relatively inexpensive AMI quantum calculation. The over- 
all expense of these more detailed charge calculations could 
be reduced by carrying them out only for the ligands with the 
best predicted binding affinity. Techniques for making con- 
strained moves of a molecule made up of rigid fragments and 
for efficiently evaluating the resulting energy (including the 
solvation free energy) are given in |il7]. 

We can further improve the accuracy of the approximate 
energy evaluation by including some contributions to SE. For 
example, we can set SE{ABC, P) « 5E{ABC) and evaluate 
SE{ABC) using eq. ( fTal together with a direct evaluation of 
E{ABC) (which, typically, is fast because it does not involve 
the interaction with the large protein). In this way, we expect 
to include the main intra-molecular contributions to the en- 
ergy including the effects of atom removal and charge redis- 
tribution. In this approximation, we still neglect three-body 
effects which enter into the solvation energy for the bound 
system, e.g., the modification of solvation energy of A inter- 
acting with P due to the presence of C |4]. Also neglected are 
the effects of atom removal and charge redistribution on the 
ligand-protein energy. 

We might make a further simplification to 5E{ABC) by in- 
cluding only some terms in energies in eq. ( fTal . For example, 
we might include just the torsion energy of the inter- fragment 
bonds and a "steric" energy, which is infinite if non-bonded 
atoms in ABC overlap and is zero otherwise. 

Finally, note that we do not need to include the chemical 
bond energies when forming ABC because these energies are 
the same in the bound and unbound systems, and so cancel in 
the computation of the binding affinity. 

C. Combining molecules 

When forming a complex ligand ABC from ligands AB 
and BC, we need to generate starting configurations for ABC 
for the purposes of identifying the wormhole portals We 



consider the problem for the bound (A = 1) case; the unbound 
case (A — 0) follows as a straightforward simplification. 

We compute Gaussian mixtures for the configurations of 
AB and BC bound, respectively, to the protein. We draw 
several sample configurations [rA,rB,rp] and [r^jPcPp] 
from each possible pair of Gaussian components selected 
from the two fits. Here Pm denotes the configuration of mol- 
ecule M. We form a "bond" constraint term, 

D'^\\rB-r'^f + \\rp-r'p\\\ 

where ||Pm — PmII is some suitable measure of the separation 
of the two configurations of molecule M and we take D > 0. 
If B is null, then ||Pb — Pg is replaced by a constraint term 
for the new bond between A and C, for example, an appro- 
priately weighted sum of the squared deviations of the bond 
length and bond angles from their ideal values. In forming D, 
we weight the various contributions so that D is an approxi- 
mate distance that atoms must move to satisfy D = 

We now seek nearby configurations for the two systems 
[Pa, Pb, Pp] and [Pg, Pc, Pp] for which the bond constraint 
is zero. This is accomplished by gradually decreasing a 
"target" constraint, Dt, from the initial value of D to zero. 
For a specific Dt, we randomly perturb the configurations in 
such a way as to meet the target constraint, D < Dt, and 
accept the new configuration with a Boltzmann probability 
min[l, exp(— Ai?/(/cT'))], where AE is the change in the 
(fit) energy and T' is an annealing temperature. In this way, 
we attempt to minimize the energy of the combined system 
subject to the bond constraint. This procedure is attempted 
at the physiological temperature T' — T and then at succes- 
sively higher temperatures until either we achieve 13 = or 
an upper temperature, T' — 2T, is reached. 

If the bond constraint can be satisfied and if we include 
the steric term for ABC in the approximate binding affinity 
calculation, we repeat the above procedure to satisfy a steric 
constraint S* = 0, where S measures the degree of overlap be- 
tween the non-bonded atoms of ABC. In this case, we perturb 
the molecule ABC subject to the constraint D = in order 
to reduce S to zero following a similar strategy as that used to 
meet the bond constraint. 

If the bond and (if applicable) steric constraints can be sat- 
isfied in such a way that both [Pa, Pb, Tp] and [Pb, Pc, Tp] 
are within a few standard deviations of one of the components 
of their respective Gaussian mixtures, then [Pa, Pb, Pc, Tp] 
gives a configuration for ABC bound to P which serves as one 
of the starting points for finding portals for the bound system. 

The scheme described above is appropriate when we have 
direct calculations of the fit energies of AB and BC interacting 
with the protein. However, if these are given by summing the 
contributions over fragments, then the bond minimization of 
the new bonds is carried out allowing all the inter-fragment 
bonds to relax. In this case, the "old" inter-fragment bonds 
would start in the ideal constrained state; however, in allowing 
the new bonds to relax, the old bonds are allowed to stretch 
so that the ligand can find a good pose where all the inter- 
fragment bonds meet the constraints. 



10 



D. Implementation 

The preceding sections describe the physical basis for using 
binding free energy to direct ligand design. The implementa- 
tion entails additional challenges. There is a need for book- 
keeping to associate a molecule with the smaller molecules 
out of which it was created. Care must be taken to match up 
the configurations of the molecules in order to implement the 
approximate energy evaluations. In order to avoid building the 
same molecule multiple times (e.g., putting the fragments to- 
gether in a different order), we use the USMILES representa- 
tion 1 18] as a unique tag for the molecule. (Unfortunately, this 
representation has shortcomings for our purposes. It does not, 
in fact, provide a unique representation of a molecule; e.g., 
C1C2CC2CC3CC13 and C1C2CC3CC3CC12 are two dif- 
ferent USMILES representations of the same molecule. Fur- 
thermore, USMILES only deals with the 2-dimensional struc- 
ture of a molecule and for the purposes of binding affinity, 
stereoisomers should be treated as distinct.) 

In our current implementation, we build molecules by 
adding one fragment at a time. The approximate binding affin- 
ity of the molecules is computed by summing the fit energy of 
the individual fragments and including the steric energy. Any 
intermediate molecule meeting a threshold binding affinity is 
recorded and is used as a base from which to build larger mol- 
ecules (up to a given size) in a depth-first fashion. 

We have tested this procedure by building ligands which 
bind to botulinum neurotoxin type B lll9ll starting with 35 or- 
ganic fragments. A subsequent evaluation of the binding affin- 
ity using the full energy shows that good agreement with the 
approximate binding affinity in the case of 2-fragment ligands, 
with 90% of the pairs having an approximate binding affinity 
within 1-2 log units of the full binding affinity. However the 
agreement is poor for ligands made up of 3-5 fragments. We 
attribute this to basing the approximate energy of A^-fragment 
ligands purely on the energies of single fragments. The errors 
in the use of the approximate energies may become excessive 
so that the approximate binding affinity is no longer close to 
the full binding affinity. Alternatively, it's possible that the 
approximate energy is still reasonably accurate but that the 
binding mode is slightly wrong so that using the approximate 
distributions to provide the initial portals for the full binding 
affinity calculation may be inadequate; if the binding mode is 
quite tight, the full binding affinity calculation may never find 
it. Both of these problems would be largely overcome by bas- 
ing the approximate energy for A^-fragment ligands on the fit 
energy for (A^ — l)-fragment ligands. 

V. DISCUSSION 

We have described a way to determine the approximate 
binding affinity of a ligands based on knowledge of the bind- 
ing affinity of simpler ligands and the associated equilibrium 
distributions. This procedure correctly accounts for the loss 
of entropy associated with connecting molecules together to 
form a larger molecule and allows the binding free energy to 
be used to direct the design of ligands. 



The predictive capability of our initial implementation (de- 
scribed above) is limited to ligands made up of just 2 frag- 
ments. However, we believe that this limitation could be 
removed by computing the full binding affinity of the best 
(A^ — 1) -fragment ligands and using this as the basis of 
building A^-fragment ligands. Such a scheme would allow a 
breadth-first search which would allow the search to be di- 
rected toward the molecules with the greatest binding affinity. 

In pruning intermediate molecules, we may wish to retain 
(A^ + l)-fragment molecules which have a worse binding 
affinity than their A^-fragment parents, with the expectation 
that this would enable us to build better (A^ + 2)-fragment 
molecules. This would allow non-functional linkers (with no 
intrinsic propensity to bind to the protein) to bridge between 
high-affinity functional groups. 

The technique described here covers joining fragments in 
a simply connected fashion where each added fragment at- 
taches at one point. It would be possible to generalize the 
method to allow the formation of rings. For example we might 
overlap the molecules ABC, BCD, CDA, DAB to create a 4- 
fragmentring (ABCD). 

The wormhole technique and the ability to fit distributions 
of molecular configurations with Gaussian mixtures intrin- 
sically depends on the system having a "small" number of 
degrees of freedom, because we require that phase space be 
spanned by a reasonable number of samples. This limits the 
degree of flexibility that can be allowed for the protein and 
dictates the use of an implicit solvent model. On the other 
hand, because we are just interested in describing where sam- 
ples are concentrated, one might expect the method to con- 
tinue to function well as the number of degrees of freedom is 
increased to, say, 20. 

Our ability to obtain realistic results is also limited by the 
accuracy of the force field and the solvent model. There are 
several areas of concern. The force field, the methods for 
determining partial atomic charges, and the solvent models 
have all been developed largely independently, and it's not 
clear how consistent these models are. It is also noteworthy 
that the data validating the GB models is based on compar- 
isons with solutions to the Poisson-Boltzmann equation |20|, 
which requires specification of the atomic radii, or is based 
on comparisons with experimental data for the absolute sol- 
vation energy (from vacuum to solution) 1 5 ] . More relevant 
for our purposes would be a comparison against experimental 
data for the changes in solvent free energy on molecular as- 
sociation. Another area of uncertainty is the charge state of 
the protein, which may have a small effect when differences 
in binding free energy are being computed (e.g., with the free 
energy perturbation method) but which may have a large ef- 
fect on the absolute binding free energy. The impact of salt 
effects is easily incorporated into GB models f2V\ . More in- 
teresting would be a principled treatment of the protonation of 
charged residues in the protein, proton exchange between the 
protein and ligand, and tautomerization of the protein or the 
ligand. The wormhole method offers a natural vehicle for such 
a treatment avoiding the need to treat protonation as a contin- 
uous process |22| and avoiding the need to add an uncharged 
ghost proton [23]. 
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